program DG_ode
    use Gauss_int
    implicit none
    real(kind=qp) :: a,b,h,x_left,x_right,x_mid
    integer :: k,n,cell_num,i,j
    !存储高斯积分点值
    real(kind=qp), dimension(:,:), allocatable :: b_value
    real(kind=qp), dimension(:,:), allocatable :: db_value
    !存储基函数在边界点处的值
    real(kind=qp), dimension(:,:), allocatable ::basis
    real(kind=qp), dimension(:,:), allocatable ::dbasis
    ! 计算刚度矩阵A
    real(kind=qp), dimension(:,:), allocatable :: M
    !存储高斯积分点
    real(kind=qp), dimension(4) :: points
    !存储高斯积分点处函数值
    real(kind=qp), dimension(4) :: fun_value
    real(kind=qp), dimension(4) :: tmp
    !存储数值解，真解等
    real(kind=qp), dimension(:), allocatable ::initial_value
    real(kind=qp), dimension(:), allocatable :: flux_term
    real(kind=qp), dimension(:,:), allocatable :: u_value
    real(kind=qp), dimension(:), allocatable ::num_value
    real(kind=qp), dimension(5) :: error
    !存储系数
    real(kind=qp), dimension(:,:),allocatable :: a_coef
    k = 2
    allocate(b_value(k+1,4))
    allocate(db_value(k+1,4))
    allocate(basis(k+1,3))
    allocate(dbasis(k+1,3))
    allocate(M(k+1,k+1))
    allocate(flux_term(k+1))
    ! 手动生成刚度矩阵
    select case(k)
        case(1)
            M(1,1) = 0.5_qp
            M(1,2) = -0.5_qp
            M(2,1) = 0.5_qp
            M(2,2) = 0.5_qp
        case(2)
            M(1,1) = 0.5_qp 
            M(1,2) = -0.5_qp 
            M(1,3) = 0.0_qp
            M(2,1) = 0.5_qp 
            M(2,2) = 0.0_qp 
            M(2,3) = -0.5_qp
            M(3,1) = 0.0_qp 
            M(3,2) = 0.5_qp 
            M(3,3) = 0.5_qp
        end select
    ! 取出高斯积分点，以及基函数在高斯点处的值
    call get_points(points)
    call compute_value(k, b_value, db_value)
    call compute_basis(k,basis,dbasis)
    !外层循环控制h的缩放次数
    a = 0.0_qp
    b = 1.0_qp
    do n= 1,5
        h = 1.0_qp / (8*(2**n))
        cell_num = 8*(2**n)
        allocate(initial_value(cell_num+1))
        allocate(a_coef(k+1,cell_num))
        allocate(u_value(cell_num,4))
        allocate(num_value(cell_num))
        call init_array(cell_num+1, initial_value)
        !内层循环控制单元数
        do i = 1, cell_num
            !获取当前单元的左右端点
            x_left = a + (i-1)*h
            x_right = a + i*h
            x_mid = (x_left + x_right) / 2.0_qp
            ! 计算f在高斯积分点处的值，注意投影到参考单元[-1,1]上
            do j = 1, 4
                tmp(j) = f(x_mid + 0.5_qp*h*points(j))
            end do 
            ! 计算flux_term
            do j = 1,k+1
                call vector_mul(4,b_value(j,:), tmp, fun_value)
                call compute_int(fun_value, flux_term(j))
                flux_term(j) = h/2.0_qp * flux_term(j)
                flux_term(j) = (flux_term(j)) + (initial_value(i)*basis(j,1))
            end do
            ! 求解系数
            a_coef(:,i) = matmul(M, flux_term)
            ! 更新u的值
            initial_value(i+1) = dot_product(basis(:,2),a_coef(:,i))
            ! 更新u_value
            do j = 1,4
                u_value(i,j) = u(x_mid + 0.5_qp*h*points(j))
            end do
        end do
        ! 计算误差
        call compute_error(h,k,cell_num,a_coef,u_value,error(n))
        !print*, error(n)
        deallocate(initial_value)
        deallocate(a_coef)
        deallocate(u_value)
        deallocate(num_value)
    end do
    call compute_order(5,error)
    !释放生成的动态数组
    deallocate(b_value)
    deallocate(db_value)
    deallocate(basis)
    deallocate(dbasis)
    deallocate(M)
    deallocate(flux_term)
    contains
    real(kind=qp) function f(x)
        implicit none
        real(kind=qp),intent(in) :: x
        f = 2*pi*cos(2*pi*x)
    end function f
    real(kind=qp) function u(x)
        implicit none
        real(kind=qp), intent(in) :: x
        u = sin(2*pi*x)
    end function
    !将数组初始化为0数组
    subroutine init_array(nn,array)
        implicit none
        integer,intent(in) :: nn
        real(kind=qp), intent(out), dimension(nn) :: array
        integer :: l
        do l= 1,nn
            array(l) = 0.0_qp
        end do   
    end subroutine init_array
end program DG_ode
include "Gauss_int.f90"